o 
o 

(N 

CO 



Universality of electron-positron distributions in extensive air showers 

S. Lafebre R. Engel ^ H. Falcke " ^ J. Horandel % T. Huege ^ J. Kuijpers \ R. Ulrich ^ 

'Department of Astrophysics, IMAPP, Radboud University, P.O. Box 9010, 6500GL Nijmegen, The Netherlands 
^Institut flir Kemphysik, Forschungszentrum Karlsruhe, P. O. Box 3640, 76021 Karlsruhe, Germany 
Radio Observatory, Astron, Dwingeloo, P.O. Box 2, 7990AA Dwingeloo, The Netherlands 



Abstract 

Using a large set of simulated extensive air showers, we investigate universality features of electron and positron distributions in very-high- 
energy cosmic-ray air showers. Most particle distributions depend only on the depth of the shower maximum and the number of particles in 
• the cascade at this depth. We provide multi-dimensional parameterizations for the electron-positron distributions in terms of particle energy, 
vertical and horizontal momentum angle, lateral distance, and time distribution of the shower front. These parameterizations can be used to 

^^btain realistic electron-positron distributions in extensive air showers for data analysis and simulations of Cherenkov radiation, fluorescence 

Q signal, and radio emission. 
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I I front structure 
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1. Introduction 

One of the greatest mysteries in particle astrophysics is the na- 
ture and origin of the highest-energy cosmic rays above 10'^ eV. 
The study of extensive air showers produced in our atmosphere 
by these particles is the primary means of obtaining information 
about high-energy cosmic rays. Many techniques to observe 
these air showers, including the detection of atmospheric fluo- 
rescence and Cherenkov light 11] and radio signal emission |2j|, 
depend on the knowledge of the distribution of charged parti- 
cles in air showers. Primarily, the distributions of electrons and 
positrons as most abundant charged particles are of importance. 
Theoretical predictions of the main production and energy loss 
processes in electromagnetic showers have been available for a 
long time [3 _4J- Modern Monte Carlo techniques greatly en- 
hance the accuracy of these estimates and allow us to calculate 
the electron-positron distributions not only in electromagnetic 
showers but also showers initiated by hadrons. 

In this work, we use simulations to investigate electron- 
positron distributions in extensive air showers and their depen- 
dence on energy, species, and zenith angle of the primary particle 
and on the evolution stage of the shower. Previous studies have 
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shown that many distributions depend only on two parameters: 
the number of particles in the extensive air shower and the lon- 
gitudinal position in the shower evolution where this maximum 
occurs |l5]|6ll2l|8llll[!0l[II][Il-This concept, which is referred 
to as universality, allows us to develop parameterizations of the 
electron-positron distributions as a function of relevant quanti- 
ties such as energy, lateral distance, and momentum angles, in 
terms of only a few parameters. 



2. Method 

Electron and positron distributions in the atmosphere were 
studied through detailed Monte Carlo simulations. Unless speci- 
fied otherwise, extensive air shower simulations were performed 
according to the specifications below. 

All simulations were carried out using the corsika code, ver- 
sion 6.5 [13 1. We used the qgsjet-II-03 model | ,14..15il to describe 
high-energy interactions and the urqmd 1.3.1 code lfT6llT7l at 
lower energies. Electromagnetic interactions were treated by 
the EGs4 code |18|. We applied a low energy cutoff of 151 keV 
and level lO""^ optimum thinning |[l9l|20]|. The U.S. Standard 
Atmosphere ET\ fT2\ was used as atmospheric model. It should 
be noted that, because simulations for our analysis were per- 
formed using only a single nuclear interaction model, the shape 
of the distributions presented may change somewhat when dif- 
ferent models such as sibyll or qgsjet-I are employed. On the 



other hand, the distributions in proton and iron showers ex- 
hibit very good universality. Hence, the overall behaviour of 
the distributions should not change significantly. 

The standard output of corsika is a list of momenta, position 
coordinates, and arrival times of those particles that cross a 
horizontal plane representing the ground detector. This output 
format is not ideally suited for universality studies. First of 
all, particle distributions need to be calculated at many depth 
layers for each individual shower. Secondly, considering inclined 
showers, different core distances in the horizontal detector plane 
correspond to different shower development stages. 

A multi-purpose interface called coast (Corsika Data Access 
Tools) has been developed for accessing the data of individual 
particles tracked in corsika fTS\. For each track segment of a 
particle simulated in corsika, a coast interface function is called 
with the particle properties at the start and end of the propagation 
step. In addition, all standard corsika output information is 
passed to the coast interface. This allows one to directly access 
the overall information of the simulated showers (e.g. energy, 
direction of incidence, depth of first interaction) as well as 
details on all individual track segments of the simulated shower 
particles. 

The COAST interface was used in this work to produce his- 
tograms of different particle distributions. Planes perpendicular 
to the shower axis were defined and particles were filled in the 
corresponding histograms if their track traversed one of these 
planes. The energy, momentum, time, and position of a particle 
crossing one of the planes was calculated by interpolation from 
the start and end points of the track segment. In total, 50 planes 
at equidistant levels in slant depth X between the point of first 
interaction and sea level (X ^ 1036 g/cm^ for vertical showers) 
were used for histogramming, whereas the depth of a plane was 
measured along the shower axis. Note that these planes are, in 
general, not horizontal and cover different atmospheric densities. 
In our universality studies below, we will use only the densities 
at the intersection points of the planes with the shower axis. 

At each of the 50 planes, two three-dimensional histograms 
were filled for electrons and positrons respectively. The first 
histogram contains logarithmically binned distributions of the 
arrival time, lateral distance from the shower axis, and the 
kinetic energy of the particles. The second histogram contains 
the angle between the momentum vector and the shower axis, 
the angle of the momentum vector projected into the plane with 
respect to the outward direction in the plane, and the kinetic 
energy of the particles. 

Showers were simulated for protons, photons, and iron nu- 
clei at primary energies of 10'^, 10'^, 10'^, and 10^" eV. For 
each combination of primary particle and energy, showers with 
zenith angles of 0, 30, 45, and 60° were calculated. Non-vertical 
showers were injected from the north, northeast, east, southeast, 
and south to accommodate deviations due to the geomagnetic 
field. The preshower effect ll24l l25l was excluded: for photon 
primaries at energies over 10^^ eV, it would result in the simula- 
tion of several lower-energy primaries, the particle distributions 
for which are already included. Each parameter set was repeated 
20 times, amounting to a total of 3840 simulated showers. The 
showers were produced with a parallelized corsika version L26il 



on a cluster of 24 nodes. Access to this library may be obtained 
through the authors. 

As a reference set, averaged distributions at the shower max- 
ima of 20 vertical air showers initiated by 10'** eV protons are 
used. This set is compared to averaged distributions of other 
parameters, only one of which is changed at a time. If not ex- 
plicitly stated, all distributions in this work refer to the sum of 
electrons and positrons. In particular, when the term 'particles' 
is used, the sum of electrons and positrons is meant. 

3. Longitudinal description 

There are several ways to describe the longitudinal evolution 
of an air shower 

Slant depth X measures the amount of matter an air shower 
has traversed in the atmosphere, in g/cm^. 

Relative evolution stage is defined here in terms of the depth 
relative to the slant depth X^ax, where the number of particles 
in the air shower reaches its maximum 



with Xq ^ 36.7 g/cm^ being the radiation length of electrons 
in air. Because the shower maximum always lies at f = 0, 
describing multiple showers in terms of this quantity rather 
than X is expected to lead to a higher degree of universality. 

Shower age is defined here so that i = at the top of the 
atmosphere, s = 1 at the shower maximum, and i = 3 at infinite 
depth 

^ _ 3X _ t + Xjnax/^O 

X + 2Xmax t/3 + Xmax/^O 

The concept of shower age arises naturally from cascade theory 
in purely electromagnetic showers |I3] |27| . For example, the 
electron energy distribution is a function of shower age. Eq. (|2|i 
is, however, only a simple, frequently used phenomenological 
approximation to the shower age parameter defined in cascade 
theory. It has the advantage that it can also be applied to showers 
with a significant hadronic component. Alternatively, shower 
age could be defined phenomenologically such that s - 
corresponds to the depth of the first interaction. Since there is 
no practical way of observing the depth of the first interaction 
in air shower measurements this variant is not considered in our 
analysis. 

To determine which description yields the highest degree of 
universality, electron energy distributions of a sample of 180 
showers of various primary energies and initiated by different 
primaries were compared. Statistical deviations from the average 
distribution were obtained at fixed relative evolution stages t 
and at each individual shower's corresponding value of X and s 
according to ([T} and ([2]). 

As an example of this comparison we show in Fig. [T] the 
statistical deviation from the mean energy distribution at each 
level. Plots are drawn as a function of f and their corresponding 
values in X and s. For descriptions in t and s, universality is 
highest near the shower maximum, because at that point all 
showers are at the same evolutionary stage by definition. This 
does not apply to the description in slant depth, where the 
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Fig. 1. Average statistical deviation from tlie average energy distribution for 
180 air sliowers of dilferent energy and primary species, averaged in slant 
depth (top), relative evolution stage (middle), and age (bottom). On average, 
the longitudinal range is the same in each plot. 

shower maxima are not lined up. In this case, the relatively fast 
evolution for younger showers is reflected in falling deviations 
with depth. When the deviation is plotted for other physical 
quantities such as momentum angle or lateral distance, all curves 
behave in a similar manner as in Fig. [T] 

Showers described in terms of X are less universal than those 
described in s or f, and slant depth is therefore rejected as 
parameter of choice. Between the two remaining descriptions, 
the diflterence is much smaller Universality is slightly better for 
descriptions in evolution stage f for f > -8, though the difference 
is insignificant. For very young showers s is a better description, 
but this stage is not of interest observationally because the 
number of particles is so small. Comparing longitudinal shower 
size profiles, if showers are compared at the same evolution 
stage f, better universality is found than when shower age s 
is used f28l. Therefore, we describe electron and positron 
distributions in terms of relative evolution stage f in this work. 

The total number of particles in the air shower crossing a 
plane at level t perpendicular to the primary's trajectory is A^(f)- 
We define 



N(t;fi) = 



dN(t) 



and n(t; fi) = 



1 dNit) 
Nif) dfi 



(3) 



as, respectively, the total and the normalised differential number 
of particles with respect to some variable fi. Likewise, distribu- 
tions as a function of two variables // and v are defined as 




Fig. 2. Average energy distribution for different evolution stages t = -6,0,6 
for electrons (marked e"), positrons (e+), and their sum (e*). Background 
curves represent simulated distributions for different primaries (p, Fe, and y) 
and energies (lO'^, lO'* and lO" eV). The corresponding parameterized 
distributions from j6} are plotted on top (dashed). 



and n(f; //, v) = 



1 d^N(t) 



(4) 



dfidv " N(t; ij.) djj. dv ' 

with dimension [yuv] ' and [v] respectively. Note that the 
definition of n(f; fj., v) implies that the distribution is normalised 
by integrating only over the last variable: 



£ 



«(f; yU, v) dv = 1, 



(5) 



making the normalisation independent of /j. In this expression, 
Vmin and v,nax are the minimum and maximum values up to 
which the histograms are calculated. 

The distributions n(t; jj, v) presented in the following sections 
may be used to obtain realistic energy-dependent particle densi- 
ties for an air shower, if the values of Xmax and A^max are given. 
One needs only to calculate the total number of particles A^(f) at 
the desired shower evolution stage. An estimate of A^(f) can be 
obtained directly from shower profile measurements or through 
one of the many parameterizations available ll29l [30l [311 l32l . 



4. Energy spectrum 

From cascade theory, the energy spectrum of electrons and 
positrons as a function of shower age takes an analytical form 



1 



as derived by Rossi & Greisen |3|; a thorough previous study of 
this parameterization was done by Nerhng et al. [T0|. Loosely 
translating this description in terms of f, we replace the equation 
by 



n(f; In e) - 



(6) 



(e + ei)r.(e + e2F^ 

where e is the energy of a given secondary particle in the shower, 
and ei,2 depend on f. We have performed a fit to this function 
for electrons, positrons and their sum, indirectly providing a 
description of the negative charge excess of extensive air showers 
as a function of evolution stage and secondary energy. In these 
fits the exponent yi was fixed at 71 -2 for positrons and 
71 = 1 for both electrons and the total number of particles. The 



parameters for all three cases are explained in Appendix A. 1 

When applied to corsika showers initiated by different species 
at different energies, the energy distribution ^ is reconstructed 
accurately. This is shown in Fig. |2] where the simulated energy 
distributions are compared to their parameterizations for evo- 
lution stages f = -6, 0, 6. For shower stages -6 < f < 9, in the 
energy region 1 MeV < e < 1 GeV, which is most relevant for 
observation of geosynchrotron or Cherenkov radiation, devia- 
tions are generally smaller than 10 % and never exceed 25 % for 
all three parameterizations. For very young showers (Fig.|2j top 
panel), increasing deviations are mainly caused by variations in 
primary energy, not by primary species type. Therefore, it high- 
lights a diminished accuracy to universally describe showers at 
f < -6 rather than hadronic model-dependence. 

Using (|6]), a similar level of universality of the energy distribu- 
tion of electrons and positrons is reached as previously obtained 
with a description in s llTOl . This basic observation is an impor- 
tant one, as it allows us to study other physical quantities in 
dependence of the electron energy in the remainder of this work. 



5. Angular spectrum 

The angular distribution of particles is an important factor 
for observations with Cherenkov and radio telescopes. For suc- 
cessful radio detection an antenna needs to be placed close to 
the shower impact position, because geosynchrotron radiation 
is beamed in a very narrow cone in the direction of propaga- 
tion |33|. As far as the particle distributions are concerned, the 
size of the patch that is illuminated on the ground then de- 
pends on the lateral distribution of the particles (cf . Sect. |7]) and 
the angle with respect to the shower axis at which they propa- 
gate. Likewise, for Cherenkov observations the angle at which 
photons are emitted is a convolution of the density-dependent 
Cherenkov angle, which is of the order of ~ 1°, and the angular 
distribution of the particles that emit them. 

Fig.[3]shows the angular distribution of particles as simulated 
in 20 individual vertical proton showers at 10'* eV as a function 
of 6. To compensate for the increase in solid angle with rising 6, 
the distribution of vertical momentum angles plotted here is 
defined in terms of Q as 
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Fig. 3. Electron distributions n{t = 0;lne, Q) at different electron energies as 
a function of momentum angle to the shower axis for 20 individual showers 
initiated by lO'* eV protons. 0° is along the primary's trajectory, 90° is 
perpendicular to the shower axis. 
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Fig. 4. Normalised average distributions n(t; In e, il) for different shower stages, 
averaged over 20 proton-initiated showers at lO"* eV. 

Since the majority of all electrons and positrons stays close to 
the shower axis, we focus on this part of the distribution. We 
will ignore the more horizontal part further away from the axis 
that can be seen at the right end of the curve for 1 GeV in Fig. [3] 
When 6 is plotted on a logarithmic scale, it becomes clear that 
there is a plateau close to the shower axis at all energies and a 
sharp drop at a certain angle that depends on secondary energy. 

Fig. |4] extends the angular distributions to different shower 
stages. The differences in the distributions are clearly smaller 
than the differences between individual showers, as noted ear- 
lier 1 10 6j JJ. The differential electron distribution with regard 
to the direction of the particle's momentum is therefore indepen- 
dent of shower stage. In addition, no perceptible dependence 
on incidence zenith angle or primary energy was found. When 
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Fig. 5. Normalised average electron distributions n{t = 0; In e, H) (solid) for 
20 proton showers at lO'* eV with 3<t statistical error margins (filled area). 
For each energy, corresponding parameterizations according to jsj are also 
drawn (dashed). 

looking at different primary species, universality seems some- 
what less convincing: spectra for heavier primary species tend 
to be wider at higher electron energies. The effect is too small, 
however, to be of consequence in our analysis. 

The universality with respect to t allows us to parameterize 
this distribution as a function of two physical quantities only: 
momentum angle and energy. We propose the form 



n(f; In e, Q.) = Co 



-l/a- 



(8) 



to describe the distribution. Values for or,- and bi, which envelop 
the dependence on e, are chosen such that the first term describes 
the flatter portion of the angular distribution parallel to the 
shower axis and the second represents the steep drop. The value 
of cr determines the smoothness of the transition from the flat 
region to the steep region. Best fit values for cr, and o-, are 



given in Appendix A.2 The dependence of these parameters 



on the secondary energy e was determined purely empirically. 
For several energies, the parameterized forms are plotted along 
with their associated simulated distributions in Fig. [5] showing 
good correspondence between the two. The parameterization 
provides a good description of the simulated distribution for the 
energy region 1 MeV < e < 10 GeV and 6 < 60°. 

We now define the cutoff angle 0c as one half of the angle at 
which e^'e"' = e''-6°'-: 



c(e) = ^ exp 



ffl — Q'2 



(9) 



For high energies, where the momentum angle is smaller than 90° 
for the majority of particles, 0^ is a measure for the root mean 
square value O^ms of the particle momentum angles. This is 
outlined in Fig. [6] in which 6^ is plotted as a function of 
energy. Theoretical root mean square scattering angles according 
to Rossi & Greisen (3) in high and low secondary energy limits 
are also drawn, as well as empirical models as parameterized 
in Hillas |,5J and Oilier et al. |,6J. At high energies, the theoretical 
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Fig. 6. Cutoft' angle 6^ according to j9j for the angular distribution as a function 
of secondary energy (solid line). Also shown are theoretical predictions 
for 

^RMs from Rossi & Greisen ||3J (dashed) as well as empirical relations 
from Hillas (5) (dash-dotted) and Giller et al. (6| (dotted). 

average scattering angle is expected oc e"', while at low energies 
it is oc e"'/^. This behaviour is reproduced properly for the cutoff 
angle. For low secondary energies (e < 3 MeV), the definition 
of a cutoff or root mean square angle becomes inapplicable 
as the angular distribution widens, covering all angles. For 
£ > 2 MeV, no appreciable difference was found between the 
angular distributions of positrons on the one hand and electrons 
on the other. 

Because our histograms do not have any sensitivity in the 
azimuthal direction by design, no dependence on the geomag- 
netic field could be determined. Previous work has shown that 
the effect on the angular distribution is probably small, but 
not negligible ||5]|34l. Because the accuracy of simulations has 
rather improved since these studies were carried out, it would 
be worthwhile to investigate the effect of the geomagnetic field 
in greater detail. 

6. Outward momentum distribution 

Let us define as the angle of a particle momentum vector 
projected in the plane perpendicular to the shower axis with 
respect to the outward direction, such that = 0° for a par- 
ticle moving away from the shower axis, and = 180° for a 
particle moving towards it. We will refer to this angle as the 
horizontal momentum angle. The effect of fluctuations in the 
horizontal angular distribution is generally much less important 
than those in the vertical angular spectrum. In fact, the distribu- 
tion of the angle of the particles does not have any influence 
on the observed signal when the distance from the observer to 
the shower is much larger than the average distance from the 
shower particles to the shower axis, as is the case in air fluores- 
cence observations. This is because the cylindrical symmetries 
of the momentum angles and the shower geometry cancel out 
independently of the shape of the distribution. Geosynchrotron 
radiation, however, will only produce a significant signal rea- 
sonably close to the shower axis, because the shower front is 
thicker in length further away (cf. Sect. |9]l, breaking down co- 
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Fig. 7. Normalised simulated horizontal angular electron distributions for 
20 individual showers initiated by lO'* eV protons at dilferent energies. 
Consecutive curve sets are shifted up by 0.005 to distinguish them better; 
curves for 1 MeV are at the actual level. 



Fig. 8. Normalised average electron distributions n(t = 0; In e, </>) (solid) for 
20 proton showers at lO'* eV with 3cr statistical error margins (filled area). 
For each energy, corresponding parameterizations according to {10) are also 
drawn (dashed). 



herence. Therefore, the horizontal momentum angle spectrum 
has to be taken into account for radio measurements. 

Simulated distributions n(f; In e, 0) at f = are plotted in 
Fig.[7]for the reference set. We observe that high-energy particles 
tend to move outward more than lower-energy particles. This 
can be explained by considering the collisions in which high- 
energy electrons and positrons are created, as they primarily 
occur close to the shower axis. Hence reaction products are 
transported away from the shower core due to their transverse 
momenta. Electrons and positrons with lower energies, on the 
other hand, are also created further away from the shower core. 

No significant dependencies on incident zenith angle, primary 
energy, and primary species were found, so the horizontal 
momentum angular spectra are universal. Additionally, the shape 
of the distribution does not change significantly for e > 2 MeV 
when only electrons or only positrons are considered. There is 
some dependence in terms of f, however: the distribution appears 
to soften with evolution stage. This effect can be explained from 
the expanding spatial structure of the shower with age. 

The distribution of n(f; 0) is very nearly exponential for 
electrons and positrons with energies over 10 GeV, while it has 
a slight bulge around the outward direction at lower energies. 
To describe the distribution, we use the parameterization 



n(t;]ne,(p) = Ci[l -i- exp(io - ^i<^ - ^24> )l 



(10) 



a form which accurately reproduces the distribution. The re- 
sulting parameter values Ao(t,e), Ai(e), and /i2(e) are explained 



in Appendix A. 3 The reference set, drawn together with its 
corresponding parameterization in Fig. [8] shows a high level of 
agreement. For other shower parameters and stages, there is a 
similar degree of consistency. 




Fig. 9. Electron distributions n{t = 0;lne, Inx) for different electron energies 
as a function of distance to the shower axis for 20 individual showers initiated 
by lO'^ eV protons. The curve set for 1 GeV is at the actual level; consecutive 
sets are shifted up by a factor of 10. 

7. Lateral distribution 

The lateral spread of particles in an air shower is of direct 
relevance since it is the primary means of obtaining informa- 
tion about the shower in ground-based scintillator experiments 
measuring particle densities at different lateral distances. By in- 
tegrating the measured distribution or using the particle density 
at a given distance, an estimate for the primary energy can be 
made. Exact knowledge of the lateral distribution shape is there- 
fore crucial to accurately determine the shape of the cosmic -ray 
energy spectrum. 

When looking at the lateral distribution of electron and 




Fig. 10. Average distributions n{t; In e, In x) for different sliower stages, averaged 
over 20 proton-initiated sliowers at lO'^ eV, clearly showing dependence on /. 
Again, consecutive sets are shifted up by a factor of 10. 

positrons in terms of tiie lateral distance r from the shower axis, 
a very poor level of universality is encountered. This is mainly 
due to differences in atmospheric density at the individual values 
of Xmax- We can compensate for these differences by expressing 
the lateral distance in terms of the Moliere unit tm, defining |35]| 



rm 9.6 g/cm^ 



(11) 



where paW is the atmospheric density as a function of height h. 
For different values of e, the normalised lateral particle distri- 
bution at / = is shown in Fig. |9] as a function of distance 
for 20 individual proton showers. In this figure, all curves line 
up as the compensation for density is applied. Note that the 
physical density A^(f ; r), expressed in particles per unit area, is 
proportional to A^(f;lnx)/x^: 



Nit;lnx) = 



dN{t) 



dA^(f) 



(12) 



d\nx ™27rrdr' 
and decreases strictly with distance from the shower axis. As 
expected, particles with higher energies tend to remain closer 
to the shower axis. This agrees with the observation that the 
angle of their momentum to the shower axis is smaller 

There is no statistically relevant dependence of the lateral 
distribution on zenith angle of incidence, nor does it change 
when electrons or positrons are considered separately, except 
at energies e < 10 MeV. There is, however, a significant effect 
with shower stage as shown in Fig. 10 older showers tend to 



be wider at the same secondary energy. Therefore, unlike in 
the case of angular distributions, in any parameterization of 
the lateral distribution a dependence on t must be incorporated. 
There is also a minor effect of the energy of the primary on the 
distribution, but this is only appreciable for secondary energies 
of e > 1 GeV. 

From Figs.|9 10 and 11 it is observed that each curve is 
a combination of two separate contributions. The left peak, 
the shape of which does not depend significantly on primary 



Fig. 1 1 . Average distributions H(f = 0; In e, In x) for different primaries, averaged 
over 20 showers at lO'* eV. Again, consecutive sets are shifted up by a factor 
of 10. Note the dependence on species of the bulge on the right. 




Fig. 12. Comparison of average distributions n{t = 0;lne, In.r) at lO'^ eV for 
20 standard photon showers to 20 proton showers in which ir- decay was 
disabled. Again, consecutive sets are shifted up by a factor of 10. 

energy or species, is produced through the main electromagnetic 
formation channel of cascading steps of bremsstrahlung and pair 
creation. The second bulge shows a high level of dependence 
on primary species, as shown in Fig. [TT[ It tends to be less 
prominent for photon primaries, as for these species there is no 
significant contribution from the pion production channel. For 
hadronic primaries it is more significant, especially at higher 
secondary energies of e > 100 MeV. The magnitude of the 
variation between different species does not change with t, but 
its lateral position does slightly. The variations in strength of 
the second bulge for different primaries can be traced back to 
the contribution initiated by the decay channel tt'^ — > //"^ -h v^. 
This is shown in Fig. 12 comparing a set of unaltered 10'^ eV 
photon-initiated showers, which have no significant pion content. 
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Fig. 13. Normalised average electron distributions n(t = 0;lne,lnx) (solid) 
for 20 proton showers at lO'* eV with 3o" statistical error margins (filled 
area). For each energy, corresponding parameterizations according to are 
also drawn (dashed). Consecutive sets are again shifted up by a factor of 10. 

to a set of proton showers at the same energy in which the 
n'^ creation channel was disabled. Differences between their 
lateral distributions are smaller than statistical deviations. 

This observation raises the question whether one could use 
this difference in lateral distribution to differentiate between pri- 
maries on an individual shower basis by their lateral distribution, 
independently of measurements of primary energy or depth of 
shower maximum. This would be a difficult task. First of all, 
appreciable difference in density only occurs at high energies 
and at some distance, implying that the total electron density 
in the region of sensitivity would be very small. Additionally, 
the effect does not appear at the same distance for different 
electron energies. This makes the feature less pronounced when 
an integrated energy spectrum is measured. 

Traditionally, the integral lateral electron distribution is de- 
scribed by a an approximation of the analytical calculation of the 
lateral distribution in electromagnetic cascades, the Nishimura- 
Kamata-Greisen (nkg) function |36, 37 1. The integral lateral 
distribution for our simulated set of showers n(f ; In x) oc jc^p^kg 
is reproduced well by a parameterization of this form, provided 
that we allow the parameters to be varied somewhat. Let us 
define 

n(f; Inx) = Cix^^xi + xY' . (13) 

as parameterization. In the original definition, described in terms 
of shower age s, we have - s, = s - 4.5, and xi - 1. 
Our simulated lateral spectra closely follow the values = 
0.0238f + 1.069, ^1 = 0.0238f - 2.918, and xi = 0.430 to an 
excellent level for 10"^ < x < 10. 

To reproduce the main bulge in the energy-dependent lateral 
electron distributions, we propose a slightly different function. 
The second bulge will be ignored here since it is much lower 
than the primary bulge, and its relative height depends heavily 
on primary species as mentioned earlier. The proposed parame- 
terization is the same as ( [T3] ): 




0.01 



e (MeV) 

Fig. 14. Cutoff distance .v^ as a function of secondary energy at different 
shower stages. The energy-independent overall break distance obtained from 
the NKG function is also plotted (horizontal line). 

«(f; In e. In x) - C2X^'«(x\ + xY' , (14) 
mimicking the behaviour of the nkg function, but now also 



varying the parameters with e. Appendix A.4 explains the values 
of x'j and As an example of the fit. Fig. 13 compares the 
parameterization to the average distribution for proton showers at 
their maximum. The proposed parameters adequately reproduce 
the main bulge of the lateral distribution in the energy range of 
1 MeV < e < 1 GeV for distances x > 2 ■ 10 and evolution 
stages -6 < f < 9. 

Neglecting the second bulge results in a slightly overestimated 
overall value for the normalisation. The disregarded tail only 
constitutes a minor fraction of the total number of particles, 
however, especially at high energies. This fact becomes even 
more evident if one considers that the actual distribution is 
obtained by dividing by x^. 

The position of the break Xc, the distance of the highest peak 
in the distribution, is plotted in Fig.[T4]for various shower stages 
for 20 averaged showers. The theoretical break distance from the 
original Nishimura-Kamata-Greisen distribution at the shower 
maximum, which is an integral distribution over all electron 
energies, is also plotted as a horizontal line. At lower energies, 
the two are in good agreement as expected. 

8. Delay time distribution 

For radio geosynchrotron measurements the arrival time of 
charged particles is a vital quantity, because it determines the 
thickness of the layer of particles that form the air shower. This 
thickness in tum defines the maximum frequency up to which 
the resulting radio signal is coherent If33ll38l . which influences 
the strength of the radio signal on the ground. 

Let us define the delay time Af of a particle as the time lag 
with respect to an imaginary particle continuing on the cosmic- 
ray primary's trajectory with the speed of light in vacuum from 
the first interaction point. In the distribution of these time lags 
we must again compensate for differences in Moliere radius to 
obtain a universal description by introducing the variable 
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Fig. 15. Election distributions n(t = 0;lne, Inr) for different electron energies 
as a function of delay time for 20 individual showers initiated by lO" eV 
protons. The curve set for 1 GeV is at the actual level; consecutive sets are 
shifted up by a factor of 10. 




Fig. 16. Average distributions n(t = 0; In 6, In t) for different primaries, averaged 
over 20 showers at lO" eV. Consecutive sets are again shifted up by a factor 
of 10. Note the species-dependent bulge on the right as in Fig. |11| 



cAt 



(15) 



where c is the speed of hght in vacuum. At sea level, t - I 
corresponds to a time delay of 0.26 fj.s. The normalised delay 
time distribution at the shower maximum for different values of e 



is shown in Fig. 15 as a function of delay time for 20 individual 



proton showers. Note the striking resemblance of the time 
lag distribution to the lateral particle distribution (cf. Fig. |9]|. 
This similarity is a direct result of the non-planar shape of the 
shower front as discussed in the next section. Therefore, every 
characteristic in the lateral distribution will have an equivalent 
in the time lag distribution. 

The dependencies on primary energy, species, and angle of 



Fig. 17. Normalised average electron distributions n(t = 0;lne, Inr) (solid) 
for 20 proton showers at lO" eV with 3o" statistical error margins (filled 
area). For each energy, corresponding parameterizations according to {16) are 
also drawn (dashed). Consecutive sets are again shifted up by a factor of 10. 

incidence closely follow those observed in the lateral distribu- 
tions in every aspect. Tliis includes the behaviour of the second 



bulge with primary species, as shown in Fig. 16 Pion-decay- 
initiated electrons and positrons are again responsible for the 
emergence of this peak. 

Given the similarity between the lateral and delay time distri- 
butions, we use a function of the same form as ( [T4] l to parame- 
terize this distribution: 



n(t; In e. In t) - C3T^»'(ti + r)^ 



(16) 



Appendix A. 5 explains the values of r, and Fig. 17 compares 



the parameterization above to the average distribution for proton 
showers at their maximum. Again, only the main peak was 
included in defining the fit parameters, causing the resulting 
parameterized shape to underestimate the number of particles 
at long delay times. 

9. Shape of the shower front 

The similarity between the lateral and delay time distributions 
of electrons and positrons as investigated in the previous sections 
is the result of the spatial extent of an air shower at a given time. 
It makes sense, therefore, to investigate the physical shape of the 
shower front by looking at the dependence of the distribution 
on lateral and delay time simultaneously. In order to keep the 
analysis practicable, we will abandon energy dependence here 
in our study. 

For 20 proton shower simulations at 10'^ eV, the shower 
front shape at the shower maximum is displayed in Fig. 18 at 



different distances from the shower core. The distribution shown 
is n(f ; In x, t), and each curve is scaled to a similar level for 
easier comparison of the distributions. Though the low number 
of particles leads to larger fluctuations of the distributions at high 
distances, the behaviour clearly does not change significantly 
for x> 3. 
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Fig. 18. Electron distributions n{t = 0;lnx,T) as a function of particle time Fig. 20. Average distributions n(t = 0;ln.v,T) for different primai'y species, 



lag for 20 individual showers initiated by lO" eV protons. 
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Fig. 19. Average distributions n{t; In x, r) for different evolution stages, averaged 
over 20 proton-initiated showers at lO" eV. 

No significant dependence of the shower front shape on 
incidence angle was found for x < 15, nor is there any change 
with primary energy. There are fluctuations with evolution 
stage, however: the time lag decreases by a constant fraction 
which depends on the shower stage. As the shower evolves, 
the entire distribution shifts to the left. This effect, shown in 



Fig. 19 can be explained from the increasing spatial structure 



of the shower with age, not unlike the case of an expanding 
spherical shell. We shall see further on that the analogy is not 
entirely legitimate, but the shift does allow one to estimate X^ax 
from the arrival times of the particles. We also found a non- 
negligible dependence of the delay time on primary species, 
which is comparable in nature to the effect of evolution stage. 



as shown in Fig. 20 The dependence of the distribution on both 
species and evolution stage can be removed almost entirely for 
distances of 0.03 < x < 15 by applying a simple exponential 



averaged over 20 proton-initiated showers at 10 eV. 

shift in T. Additionally, the distributions shown are integrated 
over energy. Therefore, the shape of the distribution changes 
when electrons or positrons are considered separately, since 
their energy distribution is different as well. 

The particle distribution at a certain distance from the shower 
core as a function of arrival time is usually parameterized as a 
gamma probability density function ll39l l40l . given by 

n(f; In X, t) oc exp[ao In T- flir]. (17) 

We have found that such a parameterization does not follow our 
simulated distributions very well. Its slope is too gentle at short 
delay times and too steep at long time lags. Here, we use the 
better representation 

n(f; Inx, t) - C4 exp[ao In t' - a\ In^ r'], (18) 

which allows for a more gradual slope on the right side of 
the curve. The modified time lag t' takes into account the 
exponential shift mentioned earlier, and is defined as 



T s re 



(19) 



where /3t and are corrections for shower evolution stage and 
primary species, respectively. The values of the parameters 
aaix), ai{x), /3t, and ySj are explained in Appendix |A.6| The 
parameter /3t can be seen as a scale width for the expansion of 
the shower front as it develops. Note that the integral lateral 
distribution as parameterized in ( [T3] l is needed to obtain actual 
particle numbers via 

A^(f ; In X, t) = N(t)n(t; In x)n(t; In x, t), (20) 

using the identities in Q. 

We may exploit the necessity of the parameter j6s in our 
description of the shower front shape to determine the primary 
species if the value of X^ax is known. To distinguish proton 
from photon showers in this manner, the required resolution in 
shower stage is 6t < y6s/y6t - 0.52, assuming perfect timing and 
distance information. This corresponds to an error in X^ax of 
19 g/cm-. To separate proton from iron showers, the maximum 
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Fig. 21. Average electron distributions n{t = 0;ln,r, r) (solid) for the refer- 
ence set with 3cr statistical error margins (filled area). For each distance, 
corresponding parameterizations according to {TsJ are drawn as well (dashed). 
Best-fit F-pdf are also plotted (dotted). 

error is reduced to 11 g/cm^. Unfortunately, these figures are 
similar to or smaller than statistical fluctuations in individual 
showers or systematic uncertainties in the atmospheric density 
due to weather influences ll4Tll42l . This makes it very difficult 
to take advantage of this intrinsic difference. 
An example of the fit of (fTS) at t 



: is shown in Fig. 21 For 
distances x > 0.8, the fit describes the simulations very accu- 
rately. Equivalence is partially lost at small distances, because 
the shape of the distribution becomes more complicated closer 
to the shower core. Even there, however, the resulting shape is 
reasonably accurate down to ^ 0.04. Also plotted are best-fit 
gamma probability density functions according to ^TT) for each 
distance, which are of lower quality than the parameterization 
used here, especially close to the core. 

For a certain distance from the shower core, we define the 
time lag as the time lag where the particle density is at its 
maximum, corresponding to the peaks of the curves shown 
earlier in this section. Its value at the shower maximum is shown 



in Fig. 22 as a function of x for the reference simulation set. 
The two straight lines represent fits of the form = Ax'' to the 
part before (dashed) and after the break (dotted) as shown in 
the plot. The time lag of the maximum particle density can be 
parameterized as 



(0.044 - 0.00170f)x' ™"° °"^^' X < xo; 
(0.028 - 0.00049f)x'-^*'"° ™""' x > xq. 



(21) 



where the value for xq follows from continuity. One could 
employ this function to estimate the value of X^^^, though the 
accuracy attainable in this way is probably much lower than 
using fluorescence measurements. 

In experiments, the shower front is sometimes approximated 
as a spherical shell |43|. How do the simulated distributions 
compare to such a hypothetical shape? Close to the shower core, 
where r <& R (with R ^ 50 the supposed curvature radius in 
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Fig. 22. Maximum density as a function of lateral distance x at the shower 
maximum. Also shown are curves for x < xq (dashed) and x > xq (dotted) 
according to the parameterization in \2l\ . 

Moliere units) we expect k - 2 and R = A Going out, the 
slope should then decrease slowly as x approaches the presumed 
curvature radius. 

This spherical shape does not correspond to the situation in 
our simulations. In the innermost region the exponent gives 
consistently smaller values of k ^ 1.79. Further out, there is 
an abrupt transition around x ^ 0.3, and the final exponent is 
1.45. 



10. Conclusion 

In this work, we have presented a framework for the accu- 
rate description of electron-positron distributions in extensive 
air showers. To characterize the longitudinal evolution of the 
air shower, the concept of slant depth relative to the shower 
maximum is used. 

Using the corsika code, we have built a library of simulations 
of air showers. Analysis of this library shows that, to a large 
extent, extensive air showers show universal behaviour at very 
high energy, making the distributions in them dependent on only 
two parameters: the atmospheric depth Xmax where the number 
of particles in the air shower peaks and the total number of 
particles A^max present in the shower at this depth. The entire 
structure of the shower follows directly from these two values. 

We have found some exceptions to the universality hypothesis 
in the spatial distribution of particles. Theoretically, these non- 
universal features can be employed to distinguish primaries on 
a shower- to- shower basis. In real experiments, however, this 
would be a difficult task because the effect either amounts to only 
a few percent, or its behaviour can be mistaken for variations 
in shower stage. 

To support the simulation of secondary radiation effects 
from extensive air showers, we have provided two-dimensional 
parameterizations to describe the electron-positron content in 
terms of stage vs. energy and stage vs. lateral distance. We have 
also supplied three-dimensional representations of the electron 
content in terms of stage vs. energy vs. vertical momentum 



angle, stage vs. energy vs. horizontal momentum angle, stage 
vs. energy vs. lateral distance close to the shower core, and 
stage vs. lateral distance vs. arrival time. 

Though these parameterizations provide accurate descriptions 
of the electron-positron distributions in air showers, the authors 
would like to mention that there are no theoretical grounds 
for most of the functional representations suggested in this 
work. Their choice is justified only by the flexibility of the 
functions to accurately reproduce the simulated distributions as 
fit functions with a small number of parameters. Additionally, 
the parameterizations provided are based on simulations with a 
single interaction model only. Though no significant changes 
are expected in the general behaviour, the parameters listed will 
likely change when a diff'erent model is employed. 

When used together with a longitudinal description for the 
total number of particles, accurate characterizations of any large 
air shower in terms of the relevant quantities can be calculated. 
These may be used for realistic electron-positron distributions 
without the need for extensive simulations and could be useful 
in calculations of fluorescence, radio or air Cherenkov signals 
from very-high-energy cosmic-ray air showers. 



A. 2. Vertical angular spectrum 

The distribution of the particles' momentum angle away from 
the shower axis can be parameterized accurately as 



n(f ; In e, Q.) = Co 



For secondary energies 1 MeV < e < 10 GeV and angles up 
to 60°, the curves are described well for n(t; In e, Q) > 10"^ by 
setting the parameters in the equations above, using nine free 
parameters, to 

bi = -3.73 + 0.92e''-2'°; 
b2 = 32.9 -4.84 In e; 

(A.l) 

ai = -0.399; 

02 = -8.36 -H 0.440 In e. 

The constant cr is a parameter describing the smoothness of 
the transition of the distribution function from the first term 
of importance near the shower axis to the second term being 
relevant further away and was set to cr = 3. The overall factor Co 
follows from the normalisation condition. 



A. 3. Horizontal angular spectrum 
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The horizontal distribution of momentum is given by 

n(t; In e,(p) = Ci[l + exp(Aa - Ai4> - A2cf>^)], ^ 

where optimal agreement is reached in the intervals 1 MeV < 
e < 10 GeV and -6 < f < 9 by setting 



^0 = 0.329 - 0.0174f + 0.669 In e - 0.0474 In^ e; 
Ai =8.10- 10"^ H- 2.79- lO^Mne; 
A2 = 1.10- lO""*- 1.14- lO^Mne, 



(A.2) 



with all energies in MeV. There were eight free parameters 
in total in the fit. The value of Ci follows directly from the 
normalisation in (|5]l. 



A.4. Lateral distribution 



Appendix A. Fit parameters 

This appendix explains in detail the various parameters used 
in the functional parameterizations throughout this paper. All 
of these were obtained by performing minimisation sequences 
using a nonlinear least-squares Marquardt-Levenberg algorithm. 



A. 1 . Energy spectrum 

The parameters in the energy spectrum distribution function 
as put forward in (|6]) were chosen to match those advocated 
in NerUng et al. [10 1. A good description is obtained with the 



parameters listed in Table A.l The constants in e\ and 62 are 
in MeV; the constant Aq is provided here for all three cases to 
obtain charge excess values; the overall parameter Ai in the 
table follows directly from normahsation constraints. 



The NKG-like function to describe the primary peak in the 
lateral distribution is defined as 

n(f; In e. In x) = C'2^'o{x\ + x)^! . 

The fit was performed in the interval 1 MeV < e < 10 GeV, 
with the additional condition that x < 5x^ in order to discard the 
second, species-dependent peak. Optimal correlation is obtained 
by using the parameters 

x\ = 0.859 - 0.0461 In' e + 0.00428 \v? e; 



= 0.0263f; 
^0 = ^^+1.34 



(A.3) 



H- 0. 160 In e - 0.0404 \n^e + 0.00276 In^ e; 



=^r-4.33. 



with nine free parameters in total. The value of e is always 
expressed in MeV. Again, the value of C'^ follows directly from 
normalisation constraints and will not be discussed here. 



n 



Table A.l 

Parameter values for the energy spectrum in j6} for species of electrons, positrons, and the sum of electrons and positrons. 





Ao 








ei 




£2 




71 


72 


Electrons 


0.485A 


exp(0.183(- 


8.17r2 


10-4) 


3.22 


- 0.0068f 


106 


- l.OOr 


1 


1 + 0.0372? 


Positrons 


0.5 16A 


exp(0.201r- 


5.42(2 


lO-'') 


4.36 


- 0.0663r 


143 


-0.15r 


2 


1 + 0.0374r 


Total 


A 


exp(0.191f- 


6.91(2 


lO-'') 


5.64 


- 0.0663f 


123 


- 0.70r 


1 


1 + 0.0374r 



A. 5. Delay time distribution 

The fit function to describe the primary peak in the delay 
time distribution is identical to that of the lateral distribution, 

«(f;lne,lnT) = C3T^«'(ti +t)^", ([16| 

and was performed at energies 1 MeV < e < 10 GeV, discarding 
the second peak. Best-fit parameters are 

Ti = exp[-2.71 + 0.0823 In e-0. 114 In^e] 
^0 = 1.70 + 0.160? -0.1421ne (A.4) 
= -3.21 

with e in MeV, using 7 free parameters in total. The constant C3 
again follows from normalisation. 

A. 6. Shape of the shower front 

The shape of the shower front is parameterized as 

n(f;lnx, r) = C4 exp[ao log t' -ailog^r'], ([TSj, 
inspired by the gamma probability distribution, with 

t' = TS-I^''-'^^ ^ 
The following parameters give optimal results: 

fl() = -6.04 + 0.707 log^ X + 0.2101og^ x 

- 0.0215 log"^ X - 0.00269 log^ x; 

(A 5) 

fli = 0.855 + 0.335 log x + 0.0387 log^ x 

- 0.00662 log^ X. 

The value for ySt is fixed at /?( - 0.20, while fi^ depends on the 
primary species: 

J0S - -0.062 for iron nuclei; 

j6s = for protons; (A.6) 

0.103 for photons. 

These parameters are valid for distances of 0.4 < x < 10^ and 
10 '* < t' < 10. 
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